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ABSTRACT 


Over  the  pasc  few  decades,  many  techniques  have  been  developed  for  Che 
numerical  solution  of  integral  equations  representing  electromagnetic  scattering 
problems.  However,  a  majority  of  these  techniques  are  limited  to  electrically 
small  scatterers,  i.e.,  below  the  resonance  range.  This  is  primarily 
because  the  amount  of  CPU  computer  time  and  storage  requirements  become  pro¬ 
hibitive  for  large  body  scatterers.  Recent  work  indicates  that  a  procedure 
based  on  the  iterative  conjugate  gradient  method  can  be  incorporated  into  con¬ 
ventional  numerical  methods  in  order  to  extend  the  range  of  application  of  the 
techniques  to  larger  geometries.  In  this  paper  we  discuss  the  conjugate  gra¬ 
dient  method  and  illustrate  several  ways  in  which  it  can  be  applied  to  electro¬ 
magnetic  scattering  problems.  The  discussion  includes  mention  of  the  advantages 
of  the  method  as  compared  to  conventional  approaches  as  well  as  some  of  Its 
limitations.  In  many  practical  scattering  problems  of  interest  at  optical  wave¬ 
lengths,  the  method  can  provide  a  convenient  means  of  treating  problems  which  are 
electrically  more  than  an  order  of  magnitude  larger  than  can  be  handled  by  other 


techniques 


1.  INTRODUCTION 


The  numerical  solution  of  Integral  equations  is  an  invaluable  technique  for 
the  treatment  of  electromagnetic  scattering  problems.  Asymptotic  techniques, 
such  as  the  Geometric  Theory  of  Diffraction  (GTD)  and  its  extensions  [1],  have 
also  been  employed  for  solving  high-frequency  scattering  problems.  However, 
this  method  is  not  well-suited  for  dealing  with  scatterers  that  cannot  be  con¬ 
veniently  described  in  terms  of  a  limited  number  of  canonical  geometries  for 
which  analytical  diffraction  coefficients  are  available.  Furthermore,  the  GTD 
is  not  very  useful  for  determining  the  near-field  behavior  and  cannot  be  applied 
to  many  practical  problems  dealing  with  inhomogeneous  or  lossy  scatterers.  An 
additional  difficulty  with  this  type  of  asymptotic  approach  is  that  it  may  be 
impossible  to  estimate  the  accuracy  of  numerical  results  based  on  such  tech¬ 
niques.  Integral  equations,  on  the  other  hand,  can  be  formulated  for  scatterers 
of  arbitrary  shape  [2],  [3],  [4].  Although  few  integral  equations  can  be  solved 
exactly,  numerical  solutions  are  readily  obtained  by  a  systematic  application  of 
techniques  such  as  the  method  of  moments  (MoM)  [3]  -  [6].  However,  in  the  past, 
the  application  of  this  approach  has  been  limited  to  bodies  that  are  electri¬ 
cally  small.  The  reason  for  this  limitation  will  be  evident  from  the  discussion 
below. 

The  moment-method  procedure  Involves  replacing  the  unknown  in  the  Integral 
equation  by  N  basis  functions  and  reduces  the  problem  to  the  solution  of  an 
N-th  order  matrix  equation.  In  most  problems,  the  unknown  density  needed  foi 
accuracy  is  at  least  ten  per  linear  wavelength.  Thus,  computer  time  and 
memory  requirements  place  an  upper  limit  on  the  size  of  the  body  to  be  analyzed. 
This  is  especially  true  at  optical  wavelengths  where  the  size  of  the  scatterer 
is  often  comparable  to  or  larger  than  the  wavelength  of  the  illumination  source. 


Current  research  in  computational  electromagnetics  includes  efforts  to  develop 
more  efficient  algorithms  for  the  solution  of  integral  equations,  thereby  per¬ 
mitting  the  analysis  of  electrically  larger  and  more  complex  geometries. 

Because  of  their  potential  for  efficiency  and  low  storage  requirements,  itera¬ 
tive  methods  are  often  incorporated  into  these  algorithms  [7]. 

In  this  paper,  we  review  one  such  iterative  technique,  viz.,  the  method  of 
conjugate  gradients  [8] ,  which  has  been  found  useful  for  the  solution  of  large 
body  scattering  problems.  Although  this  technique  was  introduced  by  Hestenes 
and  Stlefel  more  than  thirty  years  ago  for  the  iterative  solution  of  matrix 
equations,  it  has  only  recently  been  applied  extensively  to  electromagnetic 
scattering  problems  [9]  -  [19].  Recent  work  has  clarified  the  advantages  and 
the  limitations  of  the  method  and  has  shown  that  the  conjugate  gradient  method 
can,  in  many  cases,  be  used  to  solve  practical  problems  which  cannot  '.sil>  je 
treated  by  the  conventional,  direct  solution  of  matrix  equations. 

The  following  section  presents  the  conjugate  gradient  algorithm  and  iden¬ 


tifies  Important  features  of  the  method.  The  remainder  of  the  article  presents 
examples  to  illustrate  a  number  of  different  ways  in  which  the  method  can  be 
applied  in  practice. 


2.  THE  CONJUGATE  GRADIENT  METHOD 


Hescenes  and  Stiefel  introduced  the  conjugate  gradient  method  in  1952  [8]. 
Since  then,  it  has  received  much  use  in  engineering  and  applied  mathematics 
[20]  -  [23].  Although  the  method  can  be  applied  to  a  quadratic  functional  on  a 
general  Hilbert  space  H,  for  our  purposes  it  is  sufficient  to  consider  specifi¬ 
cally  Euclidean  space  E^.  For  this  space,  the  inner  product  and  norm  are 
defined  in  the  usual  manner: 

<f,g>  -  g*f  (1) 

V2 

.  ||fj|-«f,f»  (2) 

The  asterisk,  in  Equation  (1)  denotes  the  transpose-conjugate  matrix. 

We  wish  to  solve  the  linear  system 

LJ  -  E  (3) 

where  E  is  a  known  N-dimensional  vector,  J  is  an  unknown  N-dimensional  vector, 
and  L  is  a  nonsingular  N  x  N  matrix.  Let  {P^ >  be  a  set  of  N  linearly  indepen¬ 
dent  vectors  in  Ejj.  The  unknown  vector  J  can  be  represented  as 

J  *  a1P1  +  ct2P2  +  ...  +  o^Pj^  (4) 

If  the  P-vectors  are  constrained  to  satisfy  the  orthogonality  requirement 

<LA  LPt,  Pj>  -  0  i  *  j  (5) 

A 

where  L  is  the  adjoint  of  L  (the  transpose-conjugate  matrix),  the  coefficients 
in  Equation  (4)  are  given  by 


<E,LPn> 


(6) 


The  approach  embodied  in  Equations  (4)  -  (6)  is  known  as  the  conjugate  direction 
method  [8].  Note  that  in  the  literature,  this  procedure  is  often  presented  for 
the  special  case  where  L  is  a  symmetric  and  positive  definite  matrix.  We  make 
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no  such  restriction,  and  as  a  result  the  formulas  presented  here  differ  somewhat 
from  those  found,  for  Instance,  in  [20]. 

The  conjugate  gradient  method  is  a  conjugate  direction  method  which  is 
augmented  with  a  recursive  procedure  for  generating  the  P-vectors  according  to 
(5).  The  P-vectors  generated  by  the  gradient  algorithm,  are,  in  fact, 
designed  to  represent  the  solution  corresponding  to  a  specific  right-hand  side  E 
in  (3).  As  a  result,  the  solutions  of  most  equations  can  be  adequately  repre¬ 
sented  by  far  less  than  the  full  set  of  N  vectors  [23].  This  is  an  attractive 
feature  of  the  conjugate  gradient  method  over  most  other  iterative  techniques, 
providing  that  one  is  dealing  with  a  single  Incident  field. 

The  algorithm  requires  the  user  to  provide  JQ,  an  initial  estimate  of  the 
unknown  J.  Normally,  for  well-conditioned  systems,  the  choice  of  JQ  has  little 
effect  on  the  number  of  iterations  required  for  a  solution,  and  a  zero  estimate 
is  often  used.  However,  it  is  important  to  note  that  the  flexibility  of 
accepting  an  initial  estimate  permits  the  user  to  terminate  the  iteration  and 
start  anew,  using  the  current  estimate  of  the  solution  as  the  initial  estimate. 
This  feature  may  be  important  when  dealing  with  large-order  equations  in  order 
to  combat  the  round-off  errors  inevitably  introduced  due  to  the  finite  word- 
length  of  a  computer. 

A  common  form  of  the  complete  algorithm  is  given  as  follows  [24]  -  [26]: 
Initial  steps 


R  -  LJ  -  E 
o  o 

p,  -  -1A1 

1  o 


(7) 

(8) 


Iterative  steps  k  ■  1,2,... 


(9) 


<E,LPk>  _  jjL^WT 

l|Lpki!2  i|Lpki|2 


Jk  "  Jk-1  +  akPk 


(10) 


Rk  "  ^k  "  E  "  Rk-1  +  akLPk 

iiL\||2 

iilVji2 


(11) 


(12) 


.  pk+1  -  *kpk  -  L\  (13> 

The  residual  norm 

I  1*1,  I  I  I  l^lr  El  i 

tteti - n®n  (l4) 

decreases  monotonically  as  the  algorithm  progresses  and  is  useful  as  an  indica¬ 
tion  of  the  average  error  in  the  solution  after  k  steps.  The  algorithm  is 
usually  terminated  when  the  residual  norm  decreases  to  some  prescribed  value. 

The  conjugate  gradient  algorithm  can  be  derived  by  considering  the  minimi¬ 
zation  of  the  quadratic  functional  [24],  [25] 


F(Jk)  -  il\ii2 


(15) 


Specifically,  the  P-functions  are  found  by  reducing  the  gradients  of  F  at  each 
estimate  of  J  to  an  orthogonal  set  satisfying  Equation  (5).  A  set  of  such  func¬ 
tions  is  said  to  be  conjugate  with  respect  to  the  operator  iAl,  hence ,  the  name 
"conjugate  gradients." 

In  theory,  the  solution  to  an  N-dimensional  system  will  be  found  in  at  most 
N  iterations.  Even  if  no  unique  solution  exists,  the  algorithm  will  converge  in 
theory  to  the  "minimum-norm”  solution  [22]. 


Based  upon  our  own  observations,  a  numerical  solution  accurate  to  several 

digits  is  found  in  N/3  iterations  or  less  for  most  well-conditioned  systems. 

The  residual  norm  given  by  (14)  would,  for  instance,  quickly  decrease  from  its 

-4 

initial  value  to  something  less  than  1  x  10  .  For  poorly  conditioned  systems, 

the  convergence  is  naturally  slower  and  is  exacerbated  by  round-off  errors  in 
any  calculation  involving  the  operator  [22],  [23].  Under  these  conditions,  it 
is  not  uncommon  for  the  residual  norm  to  remain  relatively  constant  for  many 
iterations  or  even  to  grow  under  the  influence  of  round-off  errors.  This  situ¬ 
ation  arises  occasionally  in  practice,  and  under  such  circumstances  the  algo¬ 
rithm  may  not  converge.  This  is  true  in  spite  of  the  fact  that,  theoretically, 
the  convergence  of  the  gradient  method  is  guaranteed  in  the  absence  of  any 
machine  error.  There  are  numerous  variations  on  the  conjugate  gradient  method 
as  it  is  presented  here,  some  of  which  involve  less  computation  and  may  be  freer 
from  the  round-off  problem  [23].  This  topic  is  currently  under  investigation  by 
the  authors. 

It  is  often  necessary  in  practice  to  analyze  a  single  scatterer  for  many 
different  sources  of  illumination.  Equivalently,  it  is  desirable  to  have  an 
efficient  means  for  solving  Equation  (3)  for  many  different  right-hand  sides  E. 
If  Gaussian  elimination  is  used  to  generate  the  inverse  of  a  large  matrix,  any 
number  of  right-hand  sides  can  be  treated  at  an  almost  negligible  additional 
cost.  In  theory,  the  conjugate  gradient  method  could  be  adapted  to  test  many 
sources  efficiently  by  simultaneously  expanding  each  solution  in  terms  of  the 
single  set  (PjJ.  Unfortunately,  complications  arise  which  usually  prevent  this 
from  being  practical.  As  mentioned  earlier,  the  P-vectors  generated  for  a  given 
E  are  specifically  geared  to  represent  one  solution.  For  instance,  in  a  case 
where  the  solution  exhibited  even-symmetry,  all  of  the  P-vectors  generated  to 
represent  the  solution  were  even-symmetric.  Our  experiences  indicate  that, 


except  in  special  cases,  the  entire  set  {p^},  where  the  upper  limit  of  i  could 
be  on  the  order  of  N  for  an  N  x  N  matrix,  will  be  needed  to  treat  several  right- 
hand  sides  at  one  time.  Because  the  conjugate  gradient  method  uses  a  recursive 
algorithm  to  generate  these  functions,  round-off  errors  progressively  degrade 
their  actual  orthogonality.  In  practice,  the  severity  of  the  round-off  errors 
prevents  the  successful  generation  of  the  necessary  vectors.  For  this  reason,  it 
appears  that  a  nonrecursive  algorithm  for  generating  the  P-vectors  would  be  a 
better  candidate  for  the  treatment  of  multiple  right-hand  sides.  Any  practical 
method  for  treating  large-order  systems  would  do  so  without  the  need  to  store 
the  P-vectors  in  computer  memory.  At  this  time,  no  approach  satisfying  these 
constraints  is  known  to  the  authors,  although  they  are  currently  experimenting 
with  some  algorithms  that  might  be  suitable  for  this  purpose.  Of  course,  until 
such  a  scheme  is  available,  the  conjugate  gradient  method  may  still  be  used  to 
treat  each  right-hand  side  independently. 

Certain  Integral  equations  suffer  from  "resonances'*  which  occur  when  the 
equation  permits  homogeneous  solutions,  such  as  at  frequencies  where  a  scatterer 
is  also  a  resonant  cavity  [27]  -  [29].  When  this  occurs,  the  operator  is  ill- 
conditioned,  and  the  numerical  solution  i3  corrupted  by  the  presence  of  the 
homogeneous  solution  and  by  round-off  errors  in  the  solution  process.  The  con¬ 
jugate  gradient  method  can  be  useful  as  a  flag  to  identify  this  situation,  as 
the  convergence  of  the  method  is  significantly  slower  when  the  operator  is 
poorly  conditioned  [23].  Methods  for  circumventing  the  difficulties  arising  in 
the  "resonance”  situation  have  been  discussed  elsewhere  [27]  -  [29]. 

The  above  discussion  is  a  brief  introduction  to  the  conjugate  gradient 
method.  Readers  desiring  additional  information  are  encouraged  to  consult 
references  [8],  [10],  [11],  [20]  -  [26],  In  the  following  section,  several 
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3.  APPLICATIONS 


The  conjugate  gradient  algorithm  can  be  applied  directly  to  any  discrete 
operator  equation.  The  first  step  in  a  solution  process  is  to  select  a 
discretization  scheme  and  apply  it  to  the  integral  equation.  The  method  of 
moments  is  a  general  technique  for  converting  continuous  equations  to  matrix 
equations  [5],  and  is  used  extensively  in  computational  electromagnetics. 
However,  the  moment-method  matrix  equation  for  large-body  scattering  problems 
often  exceeds  the  limits  of  computer  memory  and  iterative  methods  which  do  not 
require  the  N  x  N  matrix  to  be  stored  offer  advantages.  The  first  example  to 
follow  illustrates  the  conjugate  gradient  solution  of  a  matrix  equation  with 
emphasis  on  minimal  storage  requirements. 

Since  the  integral  equations  of  interest  are  often  convolutional  in  form, 
the  Fast  Fourier  Transform  (FFT)  algorithm  can  in  many  cases  be  adapted  to  per¬ 
form  the  convolution.  Equations  written  as  circular  (periodic)  or  linear  con¬ 
volutions  can  be  solved  exactly  using  the  FFT,  and  at  a  considerable  savings 
over  the  comparable  matrix  multiplication.  Examples  of  these  situations  are 
given  below. 

An  alternate  approach  is  to  use  the  FFT  as  an  approximation  to  the  con¬ 
tinuous  Fourier  transform,  in  order  to  compute  convolution  integrals  which  can¬ 
not  be  easily  put  into  the  form  of  discrete  convolutions.  This  procedure  is 
also  discussed  below. 

Although  these  techniques  can  be  applied  to  any  size  problem,  we  are  pri¬ 
marily  Interested  in  geometries  which  are  electrically  too  large  to  handle  with 
conventional  matrix  methods.  Comparisons  between  the  execution  times  of  the 
iterative  methods  and  conventional  matrix  methods  have  been  presented  for  small 


geometries  [15],  [18].  To  our  knowledge,  no  systematic  comparison  has  been 
published  for  large-order  equations. 

3.1.  Solution  of  Matrix  Equations 


If  the  conjugate  gradient  algorithm  is  terminated  after  several  digits  of 
accuracy  are  obtained  in  the  solution,  it  is  often  comparable  to,  and  sometimes 
better  in  efficiency  than,  the  Gaussian  elimination  solution  of  a  single  matrix 
equation  [26].  In  the  elimination  process,  the  N  x  N  matrix  is  altered  and  must 
be  stored  in  computer  memory  or  in  an  out-of-core  peripheral  unit.  This  places 
an  upper  limit  on  the  size  of  a  scatterer  which  can  be  easily  treated  by  matrix 
methods.  The  conjugate  gradient  method  can  be  used  to  extend  this  limit,  pro¬ 
vided  that  the  matrix  elements  are  sufficiently  redundant  or  can  be  generated 
from  a  relatively  simple  function  subroutine.  Each  matrix  element  is  needed 
twice  per  iteration;  therefore,  this  approach  will  only  be  practical  for  large 
systems  if  the  elements  can  be  generated  efficiently. 

Several  examples  which  were  developed  for  conventional  matrix  solution  and 
are  well-suited  for  the  Iterative  scheme  Include  two-dimensional  dielectric 
scatterers  [30],  [31],  two-dimensional  perfectly  conducting  scatterers  [32],  and 
three-dimensional  dielectric  bodies  [33].  All  of  these  are  based  on  moment- 
method  formulations  using  so-called  pulse  basis  functions  and  point-matching, 
and  all  the  matrix  elements  are  closed-form  expressions.  For  illustration,  con¬ 
sider  TM  wave  scattering  by  perfectly  conducting  two-dimensional  cylinders,  such 
as  discussed  by  Harrington  [32].  The  matrix  elements  are  given  by 
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where  k  is  2ir/X,  R  is  the  distance  between  the  center  of  the  mth  and  nth 
mn 


subdivisions  of  the  cylinder  surface,  dQ  is  the  width  of  the  n-th  subdivision, 

(2) 

and  ( • )  is  the  Hankel  function  of  the  second  kind.  In  order  to  generate  the 
matrix  elements  efficiently,  the  Hankel  functions  are  Interpolated  from  a  look¬ 
up  table.  The  size  of  the  table  is  proportional  to  the  maximum  linear  dimension 
of  the  scatterer,  yet  requires  much  less  storage  than  the  N  x  N  matrix  it 
replaces. 

Results  for  the  conducting  cylinder  example  were  obtained  for  problems 
involving  399  x  399  matrices.  Matrix  equations  of  this  order  are  beyond  the  in- 
core  storage  capabilities  of  most  modern  computers.  Figure  1  shows  the  magni¬ 
tude  of  the  surface  current  density  on  a  square  conducting  cylinder  at  three 
frequencies.  Similar  results  can  be  obtained  for  large  dielectric  cylinders 
that  may,  in  general,  be  lossy  and  inhomogeneous.  Figure  2  shows  the  polariza¬ 
tion  current  induced  inside  a  lossy  dielectric  cylinder  consisting  of  over  500 
cells  [34].  We  note  that  van  den  Berg  has  applied  the  conjugate  gradient  method 
to  solve  a  similar  problem  involving  2500  cells  in  a  dielectric  cylinder  [11], 
[12]. 

It  is  believed  that  the  iterative  method  combined  with  the  look-up  table 
approach  will  permit  many  types  of  scattering  problems  to  be  solved  for  large 
geometries.  The  authors  have  used  it  for  the  examples  discussed  above,  and 
found  the  efficiency  satisfactory  for  large-order  equations.  In  one  case,  an 
equation  of  order  1216  was  solved  on  a  Perkln-Elmer  OS-32  computer  [34]. 

However,  not  all  integral  equations  can  be  reduced  to  matrix  equations  with 
simple  elements;  the  complexity  of  the  elements  is  usually  directly  proportional 
to  the  sophistication  (and  accuracy)  of  the  discretization  scheme  used.  As  the 
complexity  of  each  matrix  element  Increases,  the  efficiency  of  the  process 
decreases.  Yet,  even  at  its  worst,  the  method  offers  an  alternative  to  Gaussian 
elimination  for  the  solution  of  large  moment-method  matrix  equations. 
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Figure  1.  Magnitude  of  current  density  induced  on  a  square  cylinder  at  three 
frequencies.  The  cylinder  model  consisted  of  339  cells. 


3.2.  Solution  of  Periodic  Problems  Using  the  FFT 


An  example  of  a  problem  well-suited  for  solution  by  the  conjugate  gra¬ 
dient  aethod  is  the  frequency-selective  surface  (FSS)  shown  in  Figure  3. 
Frequency-selective  meshes  find  many  applications  at  optical  and  Infrared  wave¬ 
lengths  as  band  pass  filters.  The  geometry  is  an  infinite-periodic  extension  of 
a  single  rectangular  cell,  and  for  plane-wave  scattering  an  integral  equation 
for  the  vector  components  of  the  current  on  the  metallic  portion  of  the  cell  can 
be  derived  using  Floquet's  theorem  [35].  The  entire  rectangular  cell  can  be 
discretized  via  the  moment  method  to  yield  a  vector  equation  containing  scalar 
components  of  the  form 
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(18) 


Equation  (18),  written  in  the  space  domain,  represents  a  circular  discrete  con¬ 


volution.  The  Kpq  of  Equation  (18)  are  actually  each  a  double-infinite  sum¬ 


mation  over  the  Floquet  modes.  Equation  (18)  can  be  written  in  the  discrete 
Fourier  transform  domain  as 


®u8  “  JuB  Ka3 


(19) 


It  may  be  undesirable  to  compute  by  applying  the  FFT  to  the  K^j ,  and  as  an 
alternative,  the  R^  can  be  approximated  by  samples  of  the  continuous  Fourier 
transform  of  K,  which  is  usually  a  simple  analytic  expression  (not  a  summation) 
[36].  This  suggests  that  a  significant  savings  in  computation  can  be  realized 
by  operating  in  the  spectral  domain.  This  entails  using  the  FFT  to  transform 
Jnm  to  Jag,  performing  the  multiplication  indicated  in  equation  (19),  and  using 
the  inverse  FFT  to  find  from  Ea^.  This  approach  is  often  used  for  problems 
Involving  periodic  geometries  [37],  [38]. 
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Figure  3*  Frequency  selective  surface. 
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This  type  of  problem  is  ideally  suited  for  solution  by  the  conjugate  gra¬ 
dient  method  because  of  the  low  storage  requirements  imposed  by  the  K  .  and  the 

Qtp 

efficiency  of  the  FFT  algorithm.  The  approach  permits  a  wide  range  of  choices 
of  basis  and  testing  functions  in  a  moment-method  formulation  as  long  as  the 
discretization  of  the  rectangular  cell  is  evenly  spaced  in  both  dimensions. 
However,  whereas  a  conventional  moment-method  approach  is  limited  by  the  size 
matrix  that  must  be  stored,  the  conjugate  gradient  method  permits  the  treatment 
of  much  larger  problems.  Recently,  the  conjugate  gradient  method  was  used  to 
analyze  an  FSS  with  over  4000  unknowns  [17]. 

3.3.  Problems  Which  Are  Linear  Discrete  Convolutions 

The  previous  example  discussed  the  use  of  the  conjugate  gradient  method  for 

periodic  geometries.  The  circular  discrete  convolutions  that  arise  can  be 
performed  efficiently  using  the  FFT.  Certain  nonperiodic  geometries  can  be 
discretized  in  such  a  way  that  the  associated  equations  become  linear  discrete 
convolutions,  and  these  can  also  be  computed  with  the  FFT.  In  the  latter  case, 
"zero-padding"  must  be  Incorporated  into  the  algorithm  [39].  Examples  of 
problems  well-suited  to  this  approach  are  planar  structures  and  two-  and  three- 
dimensional  dielectric  bodies. 

t 

As  an  illustration,  consider  TE-wave  scattering  by  an  Imperfectly  con¬ 
ducting,  or  reslstlvely  coated  planar  strip.  Following  Ray  and  Mlttra  [18],  the 
discretized  integral  equation  can  be  written  as 
N-l 

E.  ’  V.  -  l.  ,  Jn°m-n  (20> 

n-(l-N) 

for  the  choice  of  basis  functions  pictured  in  Figure  4.  Rg  is  the  resistance  of 
the  m-th  strip  as  defined  by  Senior  [40].  Since  (20)  contains  a  discrete  con¬ 
volution,  the  FFT  can  be  Incorporated  into  a  conjugate  gradient  algorithm  to 
provide  an  efficient  solution.  Ray  and  Mittra  [18]  have  shown  that,  in  addition 


to  a  savings  in  storage,  the  conjugate  gradient  method  achieves  a  savings  in 
computation  time  as  compared  to  the  Gaussian  elimination  solution  of  Equation 
(20)  [18]. 

An  example  of  the  results  for  one  strip  configuration  is  pictured  in  Figure 
5,  following  Ray  and  Mlttra  [18].  The  three-dimensional  problem  of  scattering 
from  resistive  plates  has  been  analyzed  using  a  technique  analogous  to  that  pre¬ 
sented  for  the  resistive  strip,  and  in  particular,  used  for  plate  sizes 
requiring  over  3000  unknowns  [18],  [19]. 

3.4.  Use  of  the  FFT  in  More  General  Problems 

The  above  examples  illustrate  several  ways  in  which  the  conjugate  gradient 
method  can  be  incorporated  into  the  solution  process  for  large  scattering 
problems.  In  every  case,  the  moment-method  formulation  is  adhered  to  and 
extended  to  treat  larger  geometries.  However,  the  above  examples  all  suffer 
from  some  limitation.  The  most  general  approach  is  the  matrix  solution,  yet  it 
may  only  be  practical  if  the  matrix  elements  are  relatively  simple.  The 
FFT-based  approaches  discussed  above  are  superb  for  certain  restricted 
geometries,  but  aside  from  planar  structures  and  finite  circular  cylinders,  the 
only  nonperiodic  scatterers  which  can  be  treated  are  penetrable  dielectric 
bodies.  The  limitation  is  due  to  the  discretization  process  Itself;  in  all  the 
cases  discussed  above,  the  conjugate  gradient  algorithm  was  used  to  solve  a 
numerical  system  resulting  from  an  application  of  the  method  of  moments  to  the 
original  problem.  The  moment  method  technique  is  well-understood  and  very 
systematic;  however,  other  discretization  schemes  may  be  better  suited  for  use 
with  the  FFT. 

An  alternative  approach  to  solving  scattering  problems  is  based  upon  a  pro¬ 
cedure  known  primarily  by  the  name  "Spectral  Iterative  Technique”  [17],  [18], 
[34],  [37],  [41],  [42].  To  distinguish  this  technique,  which  involves  using  the 
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Figure  5.  Magnitude  of  the  induced  current  on  a  strip  for  several  values  of 
the  constant  surface  resistance. 
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FFT  in  an  iterative  procedure,  from  the  techniques  discussed  above  using  the  FFT 
with  the  conjugate  gradient  method,  we  introduce  the  title  "Spectral  Iterative" 
(SI)  to  those  approaches  where  a  continuous  convolution  integral  is  approximated 
using  the  FFT.  In  the  previous  examples,  the  FFT  was  used  to  perform  discrete 
convolutions.  In  the  SI  approach,  the  quantities  of  Interest  are  sampled  over  a 
regular  grid  of  points.  This  sampling  scheme  is  distinctly  different  from  the 
moment-method  discretization  scheme,  and  thus  should  be  considered  a  separate 
method.  The  conjugate  gradient  approach  discussed  in  the  examples  given  above 
is  really  nothing  more  than  a  way  to  apply  the  moment  method  to  larger  problems. 
The  conjugate  gradient  iterative  procedure  may  also  be  used  as  an  SI  technique, 
in  which  the  moment  method  formulation  is  avoided  in  favor  of  a  sampling 
discretization  that  anticipates  the  use  of  the  FFT.  The  SI  approach  was  used 
in  the  past  with  an  algorithm  based  upon  an  iteration  procedure  that  did  not 
guarantee  convergence.  The  conjugate  gradient  SI  method  does  ensure  con¬ 
vergence,  and  thus  permits  the  SI  technique  to  be  applied  to  many  additional 
scattering  problems. 

For  instance,  van  den  Berg  and  De  Hoop  have  used  the  conjugate  gradient 
method  to  analyze  scattering  from  a  rough  interface  [13].  Their  approach  incor¬ 
porated  an  FFT  approximation  to  a  Fourier  transform  integral,  and  thus  is  an  SI 
technique.  The  SI  approach  offers  a  general  method  for  the  treatment  of  large- 
body  scatterers,  as  it  appears  to  have  few  fundamental  limitations  as  to  the 
types  of  problems  that  can  be  treated  with  the  FFT.  Future  work  in  this  area 
will  produce  a  systematic  way  of  using  the  method  and  analyzing  the  accuracy  of 
the  results. 
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4.  DISCUSSION 


We  have  presented  the  conjugate  gradient  algorithm  and  discussed  some  of 
its  features.  Many  of  the  remarks  on  the  attributes  and  limitations  of  the 
method  are  based  upon  our  own  observations  made  when  using  the  method  to  solve 
pertinent  integral  equations.  In  some  cases,  our  observations  may  not  concur 
with  the  recommendations  of  others;  for  instance,  it  has  been  suggested  that  the 
conjugate  gradient  method  is  useful  for  the  solution  of  very  ill-conditioned 
equations  [10]  and  that  it  is  suited  for  handling  multiple  incident  fields. 
However,  our  findings  seem  to  be.  at  odds  with  these  suggestions.  As  the  method 
has  only  recently  been  applied  to  problems  in  electromagnetic  scattering,  it  is 
expected  that  the  future  will  provide  additional  clarification. 

The  advantages  of  the  conjugate  gradient  method  for  the  solution  of  large 
scattering  problems  lie  primarily  in  the  fact  that  in  many  cases  no  alternatives 
exist  aside  from  a  brute-force  Gaussian  elimination  solution  of  the  moment- 
method  matrix  using  out-of-core  storage.  For  many  equations,  the  conjugate 
gradient  method  permits  the  solution  of  large  problems  without  a  corresponding 
need  for  large  blocks  of  computer  memory.  The  conjugate  gradient  method  may  not 
be  the  best  choice  in  a  given  situation;  it  should  be  considered  an  alternative 
which  may  or  may  not  be  useful  depending  on  the  problem.  Needless  to  say,  much 
work  remains  in  the  extension  of  this  technique  to  general  scattering  problems. 
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